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DETERMINATION  OF  THE  HEAT  DISSIPATED  FROM  A 
SPECIMEN  UNDERGOING  CYCLIC  PLASTICITY  BY  A 
HYBRID  NUMERICAL/EXPERIMENTAL  METHOD 


INTRODUCTION 

"^When  a  metallic  specimen  is  deformed  plastically,  the  bulk  of  the  irreversible  work  is  dissipated 
in  the  form  of  heat.  It  is  generally  accepted  that  the  remaining  part  of  the  input  energy  is  consumed 
by  the  change  in  the  material’s  internal  energy.  Such  internal  energy  changes  can  be  attributed  to 
phase  changes,  development  of  residual  stresses,  translation  of  dislocations;  and  the  creation  and/or 
enlargement  of  internal  surfaces  such  as  voids.  Recent  interests  in  deformation  heating  have  primarily 
been  motivated  by  metal  forming  processes  where  the  substantial  heat  generation  greatly  influences 
the  formability^etg.,  However,  the  role  of  deformation  heating  in  fatigue  has  also  attracted 

some  attention.  *>  \ 

I-.,  I 

\  The  earliest  study  on  the  energy  dissipated  by  an  oscillating  solid  was  perhaps  by  Lord  Kelvin 

[4];  who  deduced  that  “dissipation  of  energy  is  an  inevitable  result  of  every  change  of  volume.” 
Other  early  studies  include  Hopkinson  and  Williams  [5],  and  Haigh  [6].  By  using  an  extremely  sensi¬ 
tive  extensometer  to  measure  the  length  of  steel  rods  which  had  undergone  complete  and  equal  static 
load  reversals,  the  authors  in  [5]  found  that  even  in  the  “elastic”  regime,  a  minute  permanent  defor¬ 
mation  was  detectable.  From  this,  they  estimated  the  size  of  the  stress-strain  hysteresis  area  and  thus 
the  irreversible  work.  Using  temperature  measurements  of  similar  specimens  under  fatigue  load,  they 
estimated  that  the  heat  dissipated  is  approximately  80%  of  hysteresis  energy  measured  statically.  The 
cause  of  this  discrepancy  was  thought  to  be  due  to  the  fact  that  one  test  was  done  statically  while  the 
other  was  performed  dynamically. 

As  mentioned  earlier,  we  now  know  that  the  difference  between  the  heat  dissipated  and  the  hys¬ 
teresis  energy  is  a  real  physical  phenomenon  which  is  associated  with  the  change  in  energy  state  of 
the  material.  Consideration  of  the  first  law  of  thermodynamics  quickly  leads  to  the  fact  that  if  W  is 
the  irreversible  work  rate,  and  Q  the  heat  dissipation  rate,  then  for  negligible  kinetic  effects,  the  rate 
at  which  energy  is  being  accumulated  within  the  material  U  is  given  simply  by 

U  =  W  -  Q.  (1) 

It  is  clear  that  a  material  cannot  sustain  an  indefinite  accumulation  of  energy.  It  has  been  proposed 
that  t/(=  j  U  ■  dt)  may  be  view  ,  i  ,  a  “damage-energy,”  and  that  failure  would  occur  when  U 
reaches  some  critical  value  Uc  (see  '7i  "8],  and  references  within).  The  consequence  of  the  existence 
of  Uc  can  be  enormous.  If  it  can  . .  established  that  Uc  is  indeed  a  material  constant  which  is 
independent  of  loading  and  specimen  geometry,  then  it  could  in  principle  be  used  for  monitoring  the 
residual  strength  of  a  component  in  service. 

While  this  concept  is  elegantly  simple,  the  determination  of  U  (or  Uc  for  that  matter)  is  more 
complex.  In  a  comprehensive  survey  on  the  determination  of  stored  energy  of  cold  worked  metals, 
Bever  et  al.  [9]  identified  two  basic  techniques  (Single-step  and  Two-step  methods)  used  to  determine 
the  energy  stored  in  a  cold  worked  metal.  For  the  Single-step  method,  the  work  of  deformation  and 
heat  dissipated  (generally  inferred  from  temperature  measurements)  are  simultaneously  monitored,  and 
the  stored  energy  is  determined  from  the  difference  of  these  two  quantities.  Because  such  procedures 
involve  the  difference  between  two  similarly  sized  quantities,  it  was  shown  that  small  relative  errors 
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in  each  of  the  measurements  can  lead  to  large  errors  in  the  final  result.  For  the  Two-step  methods, 
the  stored  energy  is  found  by  calorimetric  means  where  the  thermal  response  of  the  cold-worked 
specimen  is  compared  to  that  of  a  virgin  specimen.  Although  such  methods  also  involve  differences, 
those  which  incorporate  a  differential  approach  can  eliminate  a  certain  degree  of  systematic  errors. 
However,  since  measurements  in  such  techniques  are  taken  after  the  deformation,  they  could  not 
account  for  the  stored  energy  that  may  be  released  immediately  after  the  deformation.  Further,  in 
applications  to  fatigue,  while  Uc  may  in  principle  be  obtained  with  this  technique,  it  could  not  be  used 
to  monitor  the  accumulation  of  U  of  a  component  in  service. 

In  this  light,  the  Single-step  approach  is  adopted  in  this  current  work,  with  the  realization  that 
precise  measurements  of  both  W  and  Q  are  necessary.  While  W  may  be  easily  and  accurately  deter¬ 
mined  from  load-displacement  data,  accurate  measurements  of  Q  are  much  more  difficult  to  achieve. 
Most  existing  experimental  methods  are  based  either  on  the  use  of  an  electrically  heated  calibration- 
specimen  for  matching  the  temperature  measurements  with  those  of  the  actual  fatigue  specimen  (e.g., 
[5], [8], [10]),  or  by  direct  calculation  of  the  heat  passing  out  of  the  specimen  from  temperature  meas¬ 
urements  either  on  or  outside  of  the  specimen  gage  section  (e.g.,  [6],  [7]  and  [11]).  However,  all 
these  methods  require  elaborate  experimental  setups  and  cannot  adequately  account  for  any  transient 
effects.  Indeed,  most  of  the  above-mentioned  works  were  careful  in  stating  that  steady  state  condi¬ 
tions  are  required. 

Rantsevich  and  Franyuk  [12]  proposed  an  experimental/analytical  solution  for  calculating  the 
thermal  energy  losses  from  a  cyclically  loaded  specimen.  It  was  shown  that  by  knowing  temperatures 
at  two  locations  on  a  specimen,  the  steady-state  diffusion  equation  may  be  solved  for  the  source  term 
by  assuming  symmetry  and  isothermal  boundary  conditions.  In  a  later  paper  [13],  Rantsevich  esta¬ 
blished  conditions  for  steady  and  quasi-steady  assumptions  so  that  certain  transient  problems  may  be 
solved  incrementally.  However,  the  restrictions  on  symmetry  and  isothermal  conditions  appear  too 
severe  and  would  not  be  applicable  in  real  tests. 

In  this  paper,  a  hybrid  numerical/experimental  technique  for  determining  the  heat  dissipated  by  a 
cyclically  loaded  specimen  is  presented.  Treated  as  an  inverse  problem,  the  one-dimensional  diffu¬ 
sion  equation  is  solved  for  the  source  term  numerically.  The  customary  boundary  conditions  used  in 
the  direct  problems  are  replaced  by  a  least  squares  criterion  on  the  numerical  temperature  profile  to 
fit  the  measured  temperature  profile  at  discrete  time  intervals.  A  stable  scheme  is  formulated  using 
the  method  of  Lagrange  multipliers  and  finite  difference  approximations.  Because  this  method  is 
independent  of  the  temperature  boundary  conditions  and  thus  the  amount  of  heat  escaped  via  the 
grips,  standard  testing  equipment  may  be  used  with  only  minimal  requirement  for  insulating  the  test 
section  of  the  specimen  to  prevent  convective  heat  loss. 

NUMERICAL  PROCEDURE 


In  general,  the  temperature  field  of  a  solid  body  under  load  is  a  function  of  the  deformation  state 
and  the  presence  of  any  internal  and  external  heat  source  and/or  sink.  The  coupling  between  tem¬ 
perature  and  deformation  is  attributed  to  the  thermoelastic  effect,  and  the  representation  of  plastic 
energy  as  internal  heat  sources.  Because  our  interest  is  in  the  temperature  variation  from  cycle  to 
cycle,  and  not  on  its  fluctuation  within  each  cycle,  the  heat  transfer  problem  can  be  greatly  simplified 
by  considering  temperatures  which  are  averaged  over  one  complete  cycle.  Since  the  specimens  under 
consideration  are  thermally  insulated  round  bars  which  may  be  considered  thermally  thin  (Biot 
numbers  <  <  1),  the  thermal  response  may  be  described  by  the  one-dimensional  diffusion  equation: 


S6  d26  q 

St  “  dx2  PC 
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where 


0  =  temperature  averaged  over  one  cycle  (K), 
a  =  { k/pC }  =  thermal  diffusivity  (m2/sec)  , 
k  =  material  heat  conductivity  (J/m-sec-K), 
p  =  density  (kg/m3)  , 

C  =  specific  heat  (J/kg-K), 
q  =  heat  generation  rate  (J/m3-sec). 


Equation  (2)  is  usually  solved  as  a  direct  problem  in  which  given  the  source  term  q  the  tem¬ 
perature  field  is  solved  subject  to  certain  initial  and  boundary  conditions.  In  the  present  proolen:, 
however,  the  objective  is  to  find  q  (=Q)  such  that  the  temperature  field  best  fit  the  measured  data. 
Because  an  analytical  form  of  q  is  not  expected,  a  numerical  approach  to  this  problem  was  adopted. 
Equation  (2)  may  be  rewritten  in  the  Crank-Nicolson  finite  difference  form  [14],  viz: 
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where  the  subscript  i  and  superscript  t  indicate  positional  and  temporal  locations  respectively.  At  is 
the  linear  mesh  size,  and  At  is  the  time-marching  interval.  Since  the  range  of  strains  considered  are 
much  lower  than  those  necessary  to  cause  necking,  uniform  plastic  deformation  is  expected.  Hence, 
the  heat  source  q  (averaged  over  At)  is  assumed  to  be  a  function  of  time  only. 

Equation  (3)  may  be  rewritten  as 

»{_!  +  ,40!  +  0!  +  1  =  Bq  +  C,,  (4) 
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We  now  need  to  define  an  objective  function  for  which  q  may  be  solved  given  some  measured 
temperature  data.  An  obvious  choice  is  a  least-squares  criterion  on  the  difference  between  measured 
and  computed  temperature  profiles.  At  time  t,  let  there  be  measured  temperatures  7j  at  nodes  j.  The 
solution  q  is  one  which  minimizes 

S  =  L  ( 7}  -  dj)2.  (5) 

That  is,  the  solution  must  satisfy 

f-  =  0,  and  >  0.  (6) 

dq  dq 2 

One  apparently  straightforward  approach  would  be  to  solve  Eqs.  (4-6)  iteratively  by  the  algo¬ 
rithm  listed  below. 


1.  Assume  an  initial  value  for  q. 

2.  Solve  for  the  temperature  field  directly. 

3.  Evaluate  S  using  Eq.  (5),  and  estimate  q  such  that  (6)  will  be  satisfied. 

4.  Repeat  steps  2  and  3  until  Eq.  (6)  is  satisfied  within  some  specified  tolerance. 

However,  this  is  relatively  inefficient.  Furthermore,  the  solution  of  Eq.  (4)  requires  the  assignment 
of  temperature  boundary  conditions.  Although  other  researchers  have  assumed  either  adiabatic  or 
isothermal  boundary  conditions  for  solving  similar  direct  problems  (e.g.,  [2], [3]),  it  was  apparent  in 
our  preliminary  tests  that  neither  assumption  is  valid.  One  alternative  would  be  to  use  the  measured 
values  as  boundary  conditions.  However,  this  would  impose  unnecessary  bias  to  the  end  points, 
implying  no  measurement  errors  at  these  locations. 

The  above  problems  may  be  overcome  by  using  the  method  of  Lagrange  multipliers.  Treating 
Eq.  (4)  as  constraints  to  the  minimization  of  Eq.  (5),  we  write  the  function  X: 

X  =  E[(7}  -  e,)2d,  +  +  A6i  +  ei  +  l  -  Bq  -  C,)],  (7) 

i 

where 

fl  when  7}  exists, 

'  ~  "10  when  7}  does  not  exist 

and 

X,  =  Lagrange  multipliers. 


Note  that  0,  is  defined  so  that  measured  data  T,  are  not  required  for  all  i.  For  the  solution  to  exist, 
we  require 
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The  expansion  of  Eqs.  (8)  results  in  the  following  linear  set  of  equations: 
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where  all  blank  spaces  of  the  above  matrix  are  zeros.  Noting  from  Eq.  (9), 

2  b\6\  +X2  =  2b\T\ 


2  6*0*  4-  X*_j  -  2.8nTn 


(10) 


and  with  the  solution  region  confined  to  the  section  bounded  by  the  two  outermost  RTDs,  so  that 

*1  =  5*  =  1  (H) 


we  have 
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Substituting  Eq.  (12)  into  Eq.  (9)  results  in 
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The  solution  for  0,  and  q  may  thus  be  obtained  by  solving  Eq.  (13)  for  x. 

A  FORTRAN  program  was  developed  for  solving  Eq.  (13)  on  a  16  bit  microcomputer.  The 
sparseness  of  the  matrix  [M]  was  not  exploited,  and  Gaussian  elimination  was  used  in  its  inversion  for 
sake  of  simplicity.  Double  precision  was  used  throughout,  and  the  conditioning  of  [M]  was  checked 
and  found  to  be  satisfactory.  A  uniform  mesh  of  16  intervals  was  used  for  the  40  mm  (1.575  in.) 
gage  section  of  the  specimen.  This  was  selected  so  that  the  measured  temperature  locations  are  coin¬ 
cident  with  the  finite  difference  nodes.  A  marching  time  step  was  set  at  0.2  times  the  interval  for 


which  experimental  data  were  available.  For  the  1  Hz  tests,  this  was  equivalent  to  a  time  step  of  0.2 
sec  for  the  first  50  seconds,  and  then  2  seconds  thereafter.  The  required  measured  values  ( 7 ,) 
between  data  points  were  obtained  by  interpolation.  Mesh  refinements  were  performed  on  several 
runs  and  the  selected  meshes  were  found  to  be  adequate. 

EXPERIMENTAL  PROCEDURE 

Testing 

The  tests  were  performed  under  fully  reversed  sinusoidal  uniaxial  loading  in  a  closed  loop  ser- 
vohydraulic  testing  machine  operated  in  strain  control.  Two  frequencies  (1  Hz,  4  Hz)  were  con¬ 
sidered.  Temperature  and  load-strain  data  were  obtained  continuously  for  the  first  fifty  cycles  and  for 
every  tenth  cycle  thereafter.  Alignment  of  the  cross  heads  was  check  before  each  test.  A  block 
schematic  of  the  testing  setup  and  associated  instrumentation  is  shown  in  Fig.  1.  A  microcomputer 
was  used  to  control  the  test  and  instrumentation. 

Specimens  were  15.88  mm  (.625  in.),  and  19.05  mm  (.75  in.)  diameter  off  the  shelf  bar  stock 
aluminum  6061  —  T6.  No  machining  or  specimen  polishing  was  performed  before  testing.  Speci¬ 
men  length  from  grip  to  grip  was  approximately  63.5  mm  (2.5  in.)  and  the  strain  was  measured  over 
a  50.8  mm  (2.0  in.)  gage  length.  Prior  to  each  test,  the  thermal  connection  of  the  temperature 
measuring  devices  was  examined  by  applying  an  external  heat  source  to  the  bar  on  one  end  and  deter¬ 
mining  the  resulting  internal  source  term  q  numerically.  Since  there  was  actually  no  internal  source, 
q  should  be  vanishing.  For  all  cases  the  internal  source  term  obtained  was  less  than  0.01  MW/m3 
(0.268  BTU/ft3-sec),  and  the  mean  standard  error  of  the  measured  temperature  from  the  numerically 
calculated  value  was  less  than  0.01  K.  This  was  considered  satisfactory  as  the  manufacturer’s  specifi¬ 
cation  on  the  RTD  accuracy  was  of  this  order. 

Temperature  Profile  Measurement 

The  temperature  profile  of  the  specimen  was  measured  by  9  platinum  resistance  temperature 
detectors  (RTDs).  The  dimensions  of  the  RTDs  were  2.3  mm  x  2  mm  and  1  mm  thick  (0.091  in  x 
0.079  in.  x  0.039  in.).  For  the  range  of  temperatures  measured  there  exists  a  linear  relation 
between  resistance  and  temperature  [15]  given  by 

R  =  0.385  T  +  100,  (15) 

where 

T  =  temperature  (°C), 

R  =  resistance  (fl). 

The  9  RTDs  were  placed  over  a  40  mm  (1 .575  in.)  length,  and  a  thin  layer  of  silicone  paste  was 
used  to  ensure  good  thermal  contact  with  the  specimen.  The  thermal  conductivity  of  the  paste  was 
extremely  high  i.e.  140  J/m-sec-K  (16  BTU-in/hr-ft2-°F)  [16].  The  position  and  the  scan  sequence 
for  the  RTDs  for  both  frequencies  are  shown  in  Figs.  2a  and  2b.  To  minimize  convection  losses,  the 
specimen  was  wrapped  with  fiberglass  insulation. 
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The  resistance  of  the  RTDs  was  determined  by  using  4-wire  ohm  measurements.  The  tempera¬ 
ture  dat?  Aere  obtained  by  an  integrated  data  acquisition  unit  (DAU  #1)  consisting  of  a  6.5  digit 
intepr?,:iig  voltmeter,  a  20  channel  FET  multiplexer  with  a  maximum  scan  rate  of  5500  channels/sec, 
a  r  ^er  with  a  resolution  of  1  n sec,  and  memory  storage.  The  integration  time  for  each  voltage  meas¬ 
urement  was  0.1  power  line  cycle  (60  Hz).  The  RTDs  were-scanned  in  sequential  order  as  shown  in 
Fig.  2a,  2b  to  minimize  bias  toward  one  end. 

Since  the  integration  time  was  shorter  than  the  completion  of  a  60  Hz  cycle,  the  measurements 
obtained  were  susceptible  to  power  line  noise.  Therefore,  a  sampling  rate  was  chosen  to  capture  the 
alternate  rise  and  fall  of  the  60  Hz  noise  and  its  effect  could  thus  be  eliminated  upon  averaging.  This 
may  be  optimized  by  setting  the  sampling  interval  for  each  channel  to  (N  4-  1/2)  •  Tp,  where  N  is  an 
integer,  and  Tp  is  the  period  of  the  power  line  noise.  With  the  available  equipment,  N  =  2  was 
chosen  to  give  a  maximum  sampling  rate  of  24  scans/sec  for  1  Hz  and  6  scans/sec  for  4  Hz. 

Measurement  of  Hysteresis  Energy 

The  hysteresis  energy  is  given  by  the  cyclic  integral  of  the  engineering  stress  and  strain.  The 
load  and  strain  data  were  obtained  by  a  separate  data  acquisition  unit  (DAU  #  2)  with  two  high  speed 
A/D  12  bit  voltmeters  and  a  high  speed  memory  storage.  The  load  range  was  set  so  that  9000  N 
(2000  lby)  corresponded  to  a  1  volt  output  signal.  The  strain  data  were  obtained  with  an  extensometer 
with  a  50.8  mm  (2.0  in.)  gage  length  and  calibrated  so  that  a  1  volt  output  signal  corresponded  to 
0.4%  strain.  Fifty  data  points  per  cycle  were  taken  and  the  voltmeters  were  triggered  by  an  external 
pulse  occurring  at  the  beginning  of  each  cycle  from  a  function  generator  driving  the  test  (see  Fig.  1). 
The  start  of  cycle  pulses  were  also  recorded  by  a  counter  contained  in  DAU  #  2.  The  data  were 
stored  in  memory  and  later  transferred  to  a  computer  for  calculation  of  the  hysteresis  energy  by 
numerical  integration. 

System  Time  Constant 

The  RTD  response  may  be  described  by  as  a  First  order  system  and  thus  exhibits  a  time  lag  to 
an  external  stimulus.  The  system  time  constant  was  determined  by  attaching  a  RTD  to  a  specimen 
similar  to  those  to  be  tested,  and  examining  its  response  to  an  applied  cycling  elastic  load.  The  ther¬ 
moelastic  response  of  the  specimen  may  be  predicted  from  Kelvin's  law  for  adiabatic  conditions 

AT  =  -AT  A o„,  (16) 


where 


A a„  =  cyclic  amplitude  of  the  first  invariant  of  stress  (N/m2), 
K  =  thermoelastic  constant  (m2/N), 

T  =  absolute  temperature  (K), 

AT  -  cyclic  temperature  amplitude  (K). 
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The  time  constant  t  may  be  determined  from  the  measured  amplitude  and  the  amplitude 
predicted  by  Kelvin’s  law 


r  =  ^  ^(Te/Tm)2  -  1 


(17) 


where 


Tm  =  measured  response  (K), 


Te  =  predicted  thermoelastic  response  (K), 

Q  =  excitation  frequency  (sec-1). 

The  time  constant  may  also  be  obtained  independently  by  examining  the  phase  lag  <j>  between  the 
applied  load  and  the  thermoelastic  response  signals  as  given  by 


T 


-tan*. 


(18) 


Using  a  lock-in-amplifier  to  measure  the  amplitude  and  the  phase  lag  of  the  thermoelastic 
response,  the  time  constant  was  determined  to  be  less  than  0.2  sec.  by  both  equations  17  and  18. 
Numerical  simulations  indicated  that  such  time  response  had  no  significant  influence  on  the  determina¬ 
tion  of  the  heat  conversion  efficiency  for  the  tests  considered. 

RESULTS  AND  DISCUSSION 

Simulated  Tests 


To  test  the  inverse  diffusion  solution  scheme,  a  direct  problem  was  first  constructed  and  solved 
using  the  Crank-Nicolson  finite  difference  scheme.  Various  fictitious,  but  representative,  time- 
dependent  heat  source  terms,  as  well  as  asymmetric  time-dependent  Dirichlet  boundary  conditions 
were  tested.  At  1  sec  intervals,  the  temperatures  at  the  9  corresponding  RTD  positions  were 
extracted  from  the  direct  solutions.  To  simulate  actual  experimental  data,  two  types  of  errors  were 
added  to  the  numerical  outputs.  Systematic  errors  were  simulated  by  a  (<  ±  0.1  K)  random  offset 
fixed  for  each  RTD  throughout  the  run  and  random  errors  associated  with  instrumentation  noise  were 
modeled  by  a  <  ±0.1  K  signal  updated  at  each  time  increment. 

The  simulated  RTD  data  were  then  used  as  input  data  in  the  inverse  program.  After  several 
runs,  it  was  found  that  the  absolute  errors  in  the  computation  of  Q  were  relatively  independent  of  Q. 
This  means  that  the  relative  accuracy  would  improve  with  increases  in  Q.  It  was  also  noticed  that 
while  the  random  noise  determined  the  size  of  the  scattering  of  Q  about  the  mean  from  increment  to 
increment,  the  random  offsets  produced  a  scatter  of  the  mean  about  zero  from  run  to  run  (given  dif¬ 
ferent  seeds  for  the  random  errors).  From  these  runs,  it  was  found  that  for  the  material  chosen,  and 
the  above  specified  temperature  measurement  errors,  Q  may  be  determined  from  the  inversion  pro¬ 
gram  to  the  order  of  0.01  MW/m3  (0.268  BTU/ft3-sec). 
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Actual  Tests 


The  results  of  the  simulated  runs  effectively  established  the  resolution  of  our  proposed  inverse 
scheme  and  helped  in  deciding  on  the  strain  ranges  suitable  in  the  actual  fatigue  tests.  Since  the  heat 
dissipated  was  expected  to  be  near  the  hysteresis  energy,  the  strain  ranges:  ±0.522%  to  ±0.773%  at 
1  Hz  and  ±0.375%  to  ±0.571%  at  4  Hz  were  selected  to  ensure  a  hysteresis  energy  rate  of  at  least 
of  the  order  of  1  MW/m3  (26.75  BTU/ft3-sec). 

All  tests  were  run  to  complete  fracture.  Since  the  bar  stock  specimens  lacked  any  transitional 
section,  fracture  occurred  invariably  at  the  first  grip  indentation  of  one  of  the  grips.  No  apparent 
preference  to  either  the  fixed  or  moving  grip  was  noticed.  Results  for  a  typical  1  Hz  run  and  4  Hz 
run  are  shown  in  Figs.  3-12.  Figures  3  and  4  show  the  development  of  both  measured  and  computed 
temperatures  at  the  center  and  ends.  Temperature  profiles  at  three  different  times  are  also  shown  in 
Figs.  5  and  6.  It  can  be  seen  that  excellent  agreements  between  the  experimental  and  computed  tem¬ 
perature  profiles  were  achieved.  The  mean  standard  error  (MSE)  was  calculated  at  each  increment, 
and  for  all  the  specimens  tested,  the  maximum  MSE  ranged  from  0.02  K  to  0.16  K.  Such  excep¬ 
tional  agreement  validates  the  assumption  of  uniform  heat  dissipation  along  the  length  of  the  speci¬ 
men,  and  reinforces  confidence  in  its  solution. 

The  general  rise  of  the  end  RTD  readings  warrants  some  comments  on  the  treatment  of  boun¬ 
dary  conditions  by  other  analytical  and  numerical  models.  When  solving  the  direct  problem,  most 
assume  isothermal  boundary  conditions  based  on  the  argument  that  the  grips  act  as  infinite  heat  sinks. 
However,  it  is  clear  from  our  measurements  (extrapolating  to  the  grips  located  at  approximately  10 
mm  from  the  end  RTD’s)  that  such  assumptions  are  not  valid.  No  such  assumptions  were  made  in 
the  present  analysis. 

Another  interesting  feature  of  the  temperature  histories  is  the  rapid  rise  of  temperature  at  one 
end  relative  to  the  other  just  prior  to  failure.  The  fact  that  this  hotter  end  was  found  always  to  be  the 
fractured  end  suggests  that  this  rapid  increase  must  be  associated  with  the  initiation  and  subsequent 
growth  of  the  crack  since  the  intense  crack-tip  plastic  deformation  acted  as  a  large  external  (i.e.,  out¬ 
side  the  gage  length)  heat  source.  It  is  interesting  to  note  that  in  a  series  of  fatigue  tests  on  steel, 
Rantsevich  [13]  showed  that  by  analyzing  temperature  data  near  the  “vulnerable  zone,’’  the  detection 
of  crack  initiation  can  in  fact  be  made  much  earlier  than  the  magnetic  method  used  in  his  tests. 

Figures  7  and  8  show  the  comparison  of  the  computed  values  of  W  and  Q.  For  all  cases  stu¬ 
died,  the  energy  rates  showed  three  distinctive  regimes:  a  brief  but  rapid  developing  segment  at  the 
start,  followed  by  a  steady  intermediate  region  and  finally  another  rapidly  changing  section  prior  to 
complete  fracture.  The  appearance  of  these  three  regimes  was  also  found  by  Haigh  [6],  who  referred 
to  them  as  the  primary,  secondary  and  tertiary  stages,  and  attributed  the  behavior  as  a  material 
characteristic.  However,  after  careful  examination  of  the  strain  histories  from  the  current  tests,  we 
concluded  that  this  was  more  closely  associated  with  the  characteristics  of  the  testing  machine. 
Despite  testing  under  strain  control,  the  testing  machine  used  was  unable  to  maintain  a  strictly  con¬ 
stant  strain,  particular  at  the  higher  frequency  (4  Hz).  Even  for  the  1  Hz  tests  where  the  strain  his¬ 
tory  appeared  relatively  constant,  a  plot  of  the  square  of  the  strain  history  (since  near  the  elastic  limit, 
the  hysteresis  energy  is  roughly  proportional  to  the  square  of  the  strain)  revealed  the  three  distinct 
stages  resembling  those  of  the  hysteresis  history. 

The  fact  that  the  hysteresis  level  was  not  maintained  constant  throughout  the  tests  posed  no  great 
problem  to  our  present  analysis  as  steady  state  conditions  were  not  required.  In  general,  Q  was  found 
to  range  between  85%  to  95%  of  W.  Invariably,  the  heat  conversion  efficiency,  defined  as  Q/W , 
started  at  a  lower  level  and  then  gradually  rose  to  a  level  or  near-level  plateau  (see  Figs.  9  and  10). 
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At  first,  this  was  thought  to  be  due  to  an  insufficient  RTD  time  response  rate.  However,  an  analysis 
based  on  the  measured  time  constants  (see  experimental  procedure)  showed  this  effect  to  be  negligi¬ 
ble.  Another  indication  of  a  sufficient  response  time  was  the  fact  that  there  was  a  good  correlation 
between  the  fluctuations  of  W  and  Q  during  the  transient  state  (see  inset  in  Fig.  8).  Since  W  and  Q 
were  obtained  by  independent  means,  such  good  correlation  also  serves  to  support  the  correctness  in 
the  calculation  of  Q.  The  results  therefore  suggest  that  the  absorption  of  energy  (If)  is  greatest  at  the 
beginning  of  test,  and  thus  confirming  the  notion  that  the  accumulation  of  damage  is  greatest  at  the 
start  of  cycling.  For  most  cases,  Q  diverged  from  W  in  the  tertiary  stage.  The  direction  of  diver¬ 
gence  (i.e.,  whether  Q  <  W  or  W  <  Q)  appeared  unpredictable.  Since  this  tertiary  stage  invariably 
corresponded  to  the  divergence  of  the  end  RTD  readings,  indicating  crack  initiation  and  propagation, 
increasing  bending  must  have  been  present  so  that  the  extensometer  data  could  not  have  truly 
represented  the  average  strain  of  the  specimen.  This  in  turn  would  have  led  to  an  erroneous  calcula¬ 
tion  of  W.  Because  of  the  uncertainty  in  the  determination  of  IF  in  this  region,  it  was  decided  to 
exclude  it  from  the  current  analysis.  In  other  words,  we  define  failure  to  be  the  moment  at  which 
sub-critical  crack  growth  becomes  detectable  as  opposed  to  the  more  conventional  definition  of  gross 
specimen  fracture.  For  our  test  results,  this  corresponded  to  the  time  at  which  one  end  RTD  reading 
diverged  from  the  other. 

The  total  hysteresis  energy  Wc  and  the  total  absorbed  energy  Uc  up  to  the  point  of  crack  initia¬ 
tion  Nc  for  all  the  specimens  considered  are  presented  in  Fig.  11.  While  Wc  shows  a  high  depen¬ 
dency  on  Ne,  Uc  appeared  to  be  relatively  constant  over  the  cycle  range  tested.  However,  in  order  to 
confidently  establish  that  Uc  is  a  strict  constant,  more  tests  are  needed  where  a  broad  cycle  range  is 
considered  and  the  scatter  reduced.  The  use  of  higher  frequency  loads  would  permit  the  heat  dissi¬ 
pated  to  be  measured  at  lower  loads  and  thus  higher  Nc.  The  use  of  properly  designed  specimens  (as 
opposed  to  the  round-bar  stocks  used  in  the  present  work),  more  temperature  measurement  stations 
along  the  specimen,  as  well  as  a  more  precise  method  for  determining  the  moment  of  crack  initiation 
would  help  in  reducing  the  scatter  seen  in  the  present  data. 

CONCLUSION 

A  hybrid  experimental/numerical  method  for  determining  the  heat  dissipated  from  a  unifom  1/ 
deforming  specimen  is  presented.  Both  steady  and  non-steady  problems  can  be  handled  with  this  tec; 
nique,  and  it  does  not  require  any  particular  boundary  temperature  conditions  to  be  specified.  Thi 
method  was  applied  in  a  study  on  the  relationship  between  the  irreversible  work  input  and  the  heat 
dissipated  for  a  number  of  aluminum  6061-T6  specimens  under  fatigue  loading.  Its  accuracy  is  sup¬ 
ported  by  the  excellent  agreement  between  measured  and  computed  temperatures,  and  its  performance 
during  the  initial  transient  stage  is  demonstrated  by  the  good  correlation  between  the  fluctuations  in 
the  measured  hysteresis  power  W,  and  the  heat  dissipation  rate  Q  which  were  obtained  by  independent 
means.  It  was  found  that  Q  ranged  between  85%  to  95%  of  IF  and  the  rate  of  energy  absorbed  by 
the  material  U  was  invariably  greatest  at  the  start  of  the  tests,  inferring  an  initially  higher  rate  of 
damage  accumulation.  The  summing  of  energies  up  to  the  point  of  crack  initiation  revealed  that  while 
Wc  varied  greatly  with  Nc,  the  total  energy  accumulated  Uc  appeared  relatively  invariant  over  the 
range  of  cycles  considered.  However,  to  adequately  address  the  issue  on  whether  Uc  is  a  material 
parameter,  more  tests  are  needed  for  different  loading  conditions  and  higher  number  of  cycles  to 
failure. 
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Appendix  A 


PROGRAM  TSOLVE 
C 

C  program  for  solving  the  inverse  1-D  diffusion  equation.  Given 

C  temperature  readings  at  at  least  three  distinct  locations  on  a 

C  fatigue  specimen,  the  heat  source  is  solved  and  compared  to  the 

C  the  amount  of  irreversible  work  generated. 

C 

IMPLICIT  REAL* 8  (A-H.O-Z) 

PARAMETER (MMAX-100,  MXMAX-50,  MTMAX-10) 

DIMENSION  X(MMAX.MMAX) ,  Y(MMAX) ,  RHS(MMAX) 

DIMENSION  XINV(MMAX.MMAX) ,  AUX(MMAX.MMAX) 

DIMENSION  THETA(MXMAX) 

DIMENSION  XT(MXMAX) .JT(MXMAX) .JTX(MXMAX) 

DIMENSION  CC(MXMAX) 

DIMENSION  T(MTMAX),  Tl(MTMAX),  T2 (MTMAX) .  T3(MTMAX) 

DIMENSION  TEMPI (MTMAX) , TEMP2 (MTMAX) , TEMPl(MTMAX) 

DIMENSION  AT(MTMAX),  BT ( MTMAX ) ,  CT (MTMAX) 

COMMON  //TAU,  DELAY (MTMAX) 

CHARACTER*79  LINE 

CHARACTER* 15  INFIL1,  INFIL2 ,  INFIL3 

WRIT£(*, ' (1X.42HENTER  SOLUTION- SPECIFICATIONS  FILE  NAME:  .  $)') 
READ(* , ' (A) ' )  INFIL1 


WRITE(*  ' (IX. 42HENTER  TEMPERATURE-MEASUREMENTS  FILENAME:  ,  $)’) 
READ(* ,  (A) ' )  INFIL2 


WRITE(*, ' (1X.42HENTER  HYSTERYSIS  FILE  NAME:  ,  $)') 

READ(* , ' (A) ' )  INFIL3 


OPEN (UNIT-1 . FILE-INFIL1 , STATUS- ' OLD '  ) 

READ(1 , ' (A79) ' )  LINE 
READ(1,*)  NDT , NDX , XL , XK , RHO , CV 
WRITE!* , ' (/ , IX, A79) ' )  LINE 

WRITE(* ' (3X, 14 , 6X, 14 , 4X , 4G10 . 3) 1 )  NDT , NDX.XL.XK ,RHO , CV 
READ(1 , 1 (A79) ' )  LINE 
READ(1,*)  NRTD 

WRITE!*, ' (/, IX, A) ' )  'NUMBER  OF  TEMPERATURE  MEASUREMENT  POINTS' 
WRITE(*  ' (IX, 14) ' )  NRTD 
READ(1 ,  (A79) ' )  LINE 
READa,*)  (XT(I)  ,1-1, NRTD) 

WRITE(*. ' (/, 1X,A^9) ' )  LINE 

WRITE(* ' (1X.8G10.3)' )  (XT(1) ,1-1, NRTD) 

READ(1 ,  (A79^ ' )  LINE 
READ(1 , *)  TSTART 


WRITE!* ,  ' 
WRITE(* 


LINE 

TSTART 


READ! 1 ,  (A79) ' )  LINE 
READ(1 ,*)  TPULSE,  (DELAY(I) , 1-1 , NRTD) 
WRITER*, ' IX, A>9)' )  LINE 


WRITE!* 
WRITE(*,*) 


8G10 . 3) ' )  TPULSE,  (DELAY(I) ,1-1, NRTD) 


NT-NDX+1 
MX-2*NT -  3 
C 

C . NDT  number  of  time -marching  steps  between  hysterysis  data  interval 

C . NDXnumber  of  mesh  intervals  along  the  solution  length 

C . XLlength  of  specimen 

C . XKcoefficient  of  thermal  conduction 

C . RHOdensity  of  the  material 

C - CVspecific  heat  under  constant  volume 

C . NRTDnumber  of  rtd  stations 

C . TSTARTsolution  starting  time 

C . NTnumber  of  finite  difference  nodes 

C . MX  size  of  the  solution  matrix 

IF(MX  .GT.  MMAX  .OR.  NT  . GT.  MXMAX  .OR.  NRTD  . GT .  MTMAX)  THEN 
WRITE(*, ’ (A, A/A) ' )  '  ***  INSUFFICIENT  MEMORY  SPACE  ALLOCATED,  ', 
+  'CHECK  ARRAY  DECLARATIONS  PROGRAM  TERMINATED  !  ***' 

STOP 

ENDIF 


C . calculate  grid  positions  where  temperature  measurements  are  available 

IF(NRTD  . LT .  3  .OR.  NRTD  .GT.  MTMAX)  THEN 
WRITE(*, 10) 

10  FORMAT ( / '  INVALID  NUMBER  OF  TEMPERATURE  MEASUREMENTS, ' ,/ 

+  '  PROGRAM  TERMINATED ! ' ) 

STOP 

ENDIF 

C 

C . convert  dx,  xl  and  xt  to  metres,  and  delay(i)  to  time  units 
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c 

DX-0 . 001*DX 
XL-0 . 001*XL 
DO  20  I-l.NRTD 
XT( I )-0 . 001*XT( I) 

DELAY ( I ) -TPULS  E*DELAY ( I ) 

20  CONTINUE 
C 

C - -check  whether  the  rtds  are  located  sufficiently  close  to  the  FD  nodes 

DX-XL/FLOAT (NDX) 

DO  510  I-l.NRTD 
XJT-XT(I)/DX+1 
FRACTXJT-XJT - 1 NT (XJT+0 . 001 ) 

I F ( ABS ( FRACTXJT )  .GT.  0.001)  THEN 

WRITE (*, ' (/A, 12 , A) ' )  '  ***  WARNING  ***,  RTD - ' , I , '  DOES  NOT 
+  'LIE  EXACTLY  ON  A  GRID  POINT  !' 

Xl-( INT (XJT-0 . 5) - 1)*DX 
JT(I)-INT(XJT-0.5)+INT((XJT-Xl)/DX+0. 5) 

ELS  E 

JT(I )-INT (XJT+0 .001) 

ENDIF 

510  CONTINUE 


C 

C 

C 


42 

C 

45 


46 


48 

C 

C- 

C 


50 


C 


C 

C-- 

C 

C-- 

C-- 

C 

100 


c 


OPEN (UNIT-2 , FI LE-INFIL2 , STATUS-' OLD'  ) 

OPEN (UNIT-3 , FILE-INFIL3 , STATUS- ' OLD ' ) 

OPEN (UNIT-4 , FILE- ' TSOLVE .RES', STATUS- ' NEW ' ) 

obtain  first  three  rtd  measurements 

Tl(l)-0. 

DO  42  I-l.NRTD 
TEMPl(I)-0. 

CONTINUE 


READ (2 ,*)  T2(l) ,  (TEMP2(I),  I-l.NRTD) 

IF(TSTART  .GT.  T2(l))  THEN 
T1 (1)-T2 (1) 

DO  46  I-l.NRTD 
TEMP1(I)-TEMP2(I) 

CONTINUE 
GOTO  45 
ENDIF 

READ (  2 ,  * )  T3(1),<TEMP3(I), I-l.NRTD) 

DO  48  I-l.NRTD 

T1(I)-T1(1)+DELAY(I) 

T2( I)-T2(1)+DELAY(I) 

T3(l)-T3(l)+DELAY(l) 

CONTINUE 

read  in  time  and  hysterysis  rate 
TAU1-0 . 

IF(TSTART  .NE.  0.)  THEN 
CONTINUE 

READ(3 ,*)  TAUMID,  HYST 
TAU2-2 . *(TAUMID-TAU1) 

IF( (TSTART-TAU1)  .GT.  0.00001)  THEN 
TAU1-TAU2 
GOTO  50 
ENDIF 
CONTINUE 

CALL  GETRTD(T,T1,T2,T3, TEMPI ,TEMP2 , TEMP3 , AT, BT.CT.NRTD) 
CALL  FITPR0F(TEMP1 , XT ,NRTD , THETA , NT , DX) 

ENDIF 


DTO-O . 

TAU-TAU1 

B— 2  .  *DX*DX/XK 

ALPHA-XK/ ( RhO^CV ) 

QSUM-0 . 
wSUM-0 . 

DTHMAX-0 . 

HYSThysterysis  energy  rate  during  the  period  between  the 
taul  and  tau2 

SSUMtotal  heat  dissipated 
SUMtotal  irreversible  work 

CONTINUE 

?RRSUM-0 . 

READ( 3,*, END-1000)  TAUMID,  HYST 
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TAU2-2 . *TAUMID-TAU1 
C 

DTAU-TAU  2  - TAU 1 
DT-DTAU/FLOAT (NDT) 

IF(ABS((DT-DT0)/DT)  .GT.  0.001)  THEN 
A- -  2 . * ( 1+DX*DX/ (ALPHA*DT ) ) 

C-2 . *(1 . - DX*DX / ( ALPHA*DT ) ) 

C 

C . invert  matrix 

C 

IF(TAU1  .NE.  TSTART)  THEN 

WRITE(* , ' (/A/) 1 )  1  CHANGE  IN  DT.  RE- INVERTING  MATRIX  _ ' 

ELSE 

WRITE(* , ' (/A/) ' )  *  INVERTING  MATRIX  _ ' 

END  IF 

CALL  FILLMAT(A,B,MX,MMAX,NRTD,X,JT) 

CALL  MATINV(MMAX,MX ,X , XINV , AUX) 

END  IF 
C 

DO  700  JDT-l.NDT 
TAU-TAU+DT 
C 

C . interpolate  rtd  readings  at  tau 

CALL  GETRTD(T,T1,T2,T3, TEMPI, TEMP2 , TEMP3 , AT , BT , CT.NRTD) 

C 

C . 

C 

DO  620  1-2, NT-1 

CC(I-1)-C*THETA(I)  -THETA(I-l)  -THETA(I-t-l) 

620  CONTINUE 

C 

C . fill  rhs 

C 


C 

C 

C 

C 

C 


630 

640 


650 

C 


C 

C 

C 


680 


CALL  FILIJRHS(NT,MX.NRTD,CC,RHS,JT,T) 

--calculate  new  temperatures  and  q 

- -THETAcalculated  temperatures 

DO  640  1-1  NT-2 
THETA(I+l)-0. 

DO  630  J— 1  MX 

THETA( I+i ) -THETA( 1+1 ) +XINV ( I , J ) *RHS ( J ) 

CONTINUE 
CONTINUE 
XLAM2-0 . 

XLAMM-0 . 

QDOT-O . 

DO  650  J-l.MX 

XLAM2-XLAM2+XINV (NT- 1 , J ) *RHS ( J ) 

XLAMM-XLAMM+XINV (MX - 1 , J ) *RHS ( J ) 

QDOT  -QDOT  +XINV(MX  ,J)*RHS(J) 

CONTINUE 

Q-Q-K)D0T*DT 

THETA(l)  -T(l) -  0 . 5*XLAM2 

THETA ( NT ) -T ( NRTD  ) -0. 5*XLAMM 

calculate  degree  of  error  in  temperature  measurements 

DO  680  1-1,  NRTD 

DTH-THETA( JT ( I ) ) -T(I) 

IF(ABS (DTH)  .GT.  DTHMAX)  DTHMAX-ABS (DTH) 
ERRSUM-ERRSUM+DTH*DTH 
CONTINUE 


700  CONTINUE 

ERRSUM-ERRSUM/ ( NDT* ( NRTD -  2 ) ) 

ERRSTD-SQRT ( ERRSUM) 
qsuM-osuM+Q 
WSUM-WSUM+HYST*DTAU 
ETA-Q/(HYST*DTAU) 

USUM-WSUM-QSUM 

WRITE(*  ' (F10.2,  IX,  F10.3,  4(1X,G10.3))')  TAUMID,  ETA,  USUM, 

+  WSUM ,  DTHMAX,  ERRSTD 

WRITE (4, ' (1X.F10.2,  IX, F10 . 3 ,  4(1X,G10.3),  3X,  20(1H-))')  TAUMID, 
+  ETA,  USUM,  WSUM,  DTHMAX,  ERRSTD 
WRITE(4 , 1 (10(1X,F7.2)) ' )  TAU.  (THETA(I) ,  I-l.NT) 

WRITE (4 , ' ( 10( IX, F7 . 2 ) ) ' )  (T<I),  I-l.NRTD) 

TAU1-TAU 
DTO-DT 
GOTO  100 
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1000  CONTINUE 


10 

15 


20 

C 

C-- 

c 


25 


AVETA-100 . * ( 1 . -  USUM/WSUM) 
WRITER  (^X,A,G10 . 3 ,A) ' ) 

+WRITE(4/{/lX,A,G10.3,A)  * ) 
+Y  —  1  AVFTA  '  » ' 

WRITE(*.  ’^lX.A.GlO.S.A) 1 ) 


30 


•TOTAL  ACCUMULATED  ENERGY  -  * ,USUM, 
'AVERAGE  WORK-HEAT  CONVERSION  EFFICIENC 
•TOTAL  ACCUMULATED  ENERGY  -  ' ,USUM, 
•AVERAGE  WORK-HEAT  CONVERSION  EFFICIENC 


+  •  J/M  3 

WRITE!*, ' !/lX, A.G10 . 3 ,A) ' ) 

+Y  —  '  AVEiA  '  % ' 

WRITE!*, ' (/IX, A)’)  '***  RUN  COMPLETED  ***' 

STOP 

END 

SUBROUTINE  GETRTD !T , T1 , T2 , T3 , TEMPI , TEMP2 , TEMP3 , AT , BT , CT , NRTD ) 

Routine  for  interpolating  the  rtd  readings. 

Parabolic  interpolation  is  used  when  data  spacing  is  less  than  or  equal 
to  2  sec.  Otherwise,  linear  interpolation  is  used. 

IMPLICIT  REAL* 8  (A-H.0-Z) 

DIMENSION  T!NRTD),  Tl(NRTD) .  T2(NRTD) ,  T3(NRTD) 

DIMENSION  TEMP1!NRTD),  TEMP2(NRTD) ,  TEMP3(NRTD) 

DIMENSION  AT!NRTD),  BT!NRTD) ,  CT!NRTD) 

COMMON  //TAU,  DELAY! 1) 

LOGICAL  LINEAR 
DATA  LINEAR  /.FALSE./ 

IF!TAU  .GT.  T3!l)  .OR.  (AT!1)+BT( 1)+CT(1) )  .EQ.  0.)  THEN 
IF!TAU  .GT.  T3(l)+l.E-6)  THEN 
DO  10  I— 1 .NRTD 
Tim-T2(I) 

T2(I)-T3(I) 

TEMPlm-TEMP2!I) 

TEMP2 ( I ) -TEMP3 ( I ) 

CONTINUE 

READ! 2,*. END-1000)  T3!l) . (TEMP3(I) , 1-1 , NRTD) 

IF(TAU  .GT.  T3(l))  GOTO  15 
DO  20  I— 1, NRTD 

T3(I)-T3!1)+DELAY(I) 

CONTINUE 

ENDIF 


-re-evaluate  polynomial  coefficients 

2.)  THEN 


IF!T3!1)-T2!1)  .GT. 

LINEAR-. TRUE. 

DO  25  I— 1 , NRTD 

BT  - 

CT  ( 

CONTINUE 
ELSE 

LINEAR-. FALSE. 

DO  30  I-l.NRTD 

AT! I)-! !T1!I) -T3(I) )*!TEMPl(I) -TEMP 

TEMP1!I)-TEMP3!I))!  '  - 

Tim-T3(I))-!Tl!I)*Tl!l)-T3!I)*T3!I))* 

T1(I) -T2(I) ) ) 

BT!I)-?TEMPl(I) -TEMP2 (I) -AT(I)*(T1!I)*T1 (I) -T2 (I)*T2 (I) ) )/ 
,T1(I) -T2(l) ) 

CT!I)-TEMP1!I)-AT!I)*T1!I)*T1!I)-BT!I)*T1!I) 

CONTINUE 

ENDIF 

ENDIF 

-interpolate  rtd  readings 


)-TEMP2!I))-!Tl(I)-T2(I))* 

!!Tl!I)*Tl!I)-T2!I)*T2!I)) 


IF!LINEAR)  THEN 
DO  40  I-l.NRTD 

T(I)-BT(I)*TAU+CT!I) 

40  CONTINUE 

ELSE 

DO  50  I-l.NRTD 

T ! I ) -AT ! I ) *TAU*TAU+BT ! I ) *TAU+CT ! I ) 

50  CONTINUE 

ENDIF 

RETURN 

1000  WRITE!*,*)  *  END  OF  FILE  IN  READING  RTD  MEASUREMENTS ’ 
WRITE!*,*)  ’  PROGRAM  TERMINATED  !' 

STOP 

END 
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SUBROUTINE  FILLMAT(A,B,M,MMAX,NRTD,X, JT) 
IMPLICIT  REAL* 8  (A-H.O-Z) 

DIMENSION  X(MMAX.MMAX) 

DIMENSION  JT(NRTD) 

C 

NT-(M+3)/2 
DO  500  I-l.M 
DO  510  J-l ,M 
X(I,J)-0. 

510  CONTINUE 
500  CONTINUE 


520 


530 


540 

C 


550 


xa,i)-A 

x  1,2  -1. 

X(l, NT-1)— 0.5 
X(1,M)-B 
DO  520  1-2, NT-3 
xa.i-D-i. 

X  I,I)-A 
X(I,I+1)-1. 
X(I,M)— B 
CONTINUE 
X(NT- 2 ,NT - 3)-l 
X(NT-2,NT-2)-A 
X(NT-2,M-1)— 0.5 
X(NT-2,M)  — B 

X(NT-1,NT-1)-A 
X(NT- 1 ,NT)— 1 . 

DO  530  I-NT.M-2 
X(I,I-l)-i. 

X  I,I)-A 
X(I.I+1)-1. 
CONTINUE 
X(M-l,M-2)-l. 
X(M-1,M-1)-A 

DO  540  I-NT-1 ,M-1 
X(M, I)— 1 . 
CONTINUE 
X(M,M)-2 


DO  550  I-2.NRTD-1 


RETURN 

END 


SUBROUTINE  MATINV (MMAX , M , A , AINV , AA) 
IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  A(MMAX.MMAX),  AINV (MMAX, MMAX) 
DIMENSION  AA(MMAX.MMAX) 


C 

C 

C 

C500 


520 

510 


PRINT*, 'A-' 

DO  500  I-l.M 

WRITE(*,  (IX, 9G10. 3) ' )  (A(I , J) , J-1,M) 
CONTINUE 
DO  510  I— 1 ,M 
DO  520  J-l.M 
AINV(I , J)-0 . 

AA(I,J)-A(I,J) 

CONTINUE 
AINV (I , I)— 1 . 

CONTINUE 


DO  530  K-2 ,M 
DO  540  I-K.M 
DO  550  KK-l.M 

AINV(I , KK)-AINV(I ,KK)-AINV(K-1,KK)*(A(I,K-1)/A(K-1,K-1)) 
550  CONTINUE 

DO  560  J-K.M 

A(I , J)— A(I ,J)-A(I,K-1)*(A(K-1,J)/A(K-1,K-1)) 

560  CONTINUE 

540  CONTINUE 
530  CONTINUE 


DO  570  KK-l.M 

AINV(M,KK)-AINV(M,KK)/A(M,M) 

DO  580  I— M- 1 ,1,-1 
DO  590  J-I+l.M 

AINV (I ,KK)-AINV(I ,KK) -A(I , J)*AINV(J ,KK) 
590  CONTINUE 
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AINV( I ,KK)-AINV( I ,KK)/A(I , I ) 

580  CONTINUE 
570  CONTINUE 
C 

C  PRINT*, ’A"- 1-' 

C  DO  600  I-l.M 

C  WRITE (* , * ( IX , 9G10 . 3 )  1  )  (AINV( I , J ) , J-l ,M) 

C600  CONTINUE 


AMAX-0 . 

DO  610  I-l.M 
DO  620  J-l.M 
ATEST-0 . 

DO  630  K-1,M 

ATEST-ATEST+AA( I , K) *AINV (K , J ) 

630  CONTINUE 

IF(I  .EQ.  J)  ATEST-ATEST - 1 . 

IF(ABS(ATEST)  .GT.  AMAX)  AMAX-ABS (ATEST) 

620  CONTINUE 
610  CONTINUE 

IF (AMAX  .GT.  l.E-10) 

+  WRITE(* , ’ (A) ' )  '***  WARNING,  MATRIX  INVERSION  MAY  BE’, 
+  '  INACCURATE  ***' 

C  PRINT*,1 I-' 

C  DO  640  I-l.M 

C  WRITE (* ,  (IX, 9G10 . 3) * )  (A(I , J) , J-l ,M) 

C640  CONTINUE 
C 

RETURN 

END 


C 

C 

C 

C 


C 


500 


SUBROUTINE  FILLRHS (NT , M , NRTD ,  C , RHS , JT , T) 

Routine  for  filling  in  the  RHS  of  the  set  of  linear  equations  to  be 
solve . 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  C(NT),  RHS(M) ,  JT(NRTD) ,  T(NRTD) 

DO  500  1-2, NT -3 
RHS (I)-C(I) 

CONTINUE 

RHS(1)-C(1)-T(1) 

RHS(NT-2)-C(NT-2) -T(NRTD) 


510 


DO  510  1-2, NRTD- 1 

RHS (NT+JT(I) -3)-2 .*T(I) 
CONTINUE 
RETURN 
END 


SUBROUTINE  FITPROF (T , XT , M , THETA , NT , DX) 
IMPLICIT  REAL* 8  (A-H.O-Z) 

DIMENSION  T(M) ,XT(M) 

DIMENSION  THETA (NT) 

DIMENSION  X(3, 3) .  XAUX(3,3),  R(3),  RAUX(3) 


C . FIT  PARABOLA 

C 

IF(M  .EQ.  3)  THEN 
DO  520  1-1,3 

X  ( 1 , 1 ) -XT ( I ) *XT ( I ) 

X(I,2 )-XT(I) 

X(1 , 3)-l . 

R(I)-T(I) 

520  CONTINUE 

ELSE 

DO  530  1-1 ,M 
XSUM-XSUM+XT ( I ) 

X2SUM-X2SUM+XT ( I ) *XT ( I ) 
X3SUM-X3SUM+XT ( I ) *XT ( I ) *XT ( I ) 
X4SUM-X4SUM+XT ( I ) *XT ( I ) *XT ( I ) *XT ( I ) 
XTSUM-XTSUM+XT ( I ) *T ( I ) 
X2TSUM-X2TSUM+XT ( I ) *XT ( I ) *T ( I ) 
TSUM-TSUM+T ( I ) 

530  CONTINUE 

X(1 , 1)-X4SUM 
X(1,2)-X3SUM 
X(1 , 3 J-X2SUM 
X(2,1)-X3SUM 
X( 2 , 2)-X2SUM 
X(2 , 3)-XSUM 
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X(3,1)-X2SUM 
X(3,2)-XSUM 
X(3 , 3)-M 
R(1)-X2TSUM 
R(2)-XTSUM 
R(3)-TSUM 
ENDIF 

CALL  MATSOL( 3 , 3 ,X, R tXAUX , RAUX, 1) 

DO  540  I-l.NT 
XX-(I-1)*DX 

THETA ( I ) -R  < 1 ) *XX*XX+R ( 2 ) *XX+R ( 3 ) 

540  CONTINUE 
C 

RETURN 

END 

SUBROUTINE  MATSOL(M , MMAX , A , B , AA , BB , IOPT) 

IMPLICIT  REAL* 8  (A-H.O-Z) 

DIMENSION  A(MMAX.MMAX) ,  B(M) 

DIMENSION  AA  (MMAX, MMAX )  ,  BB(M) 

DO  500  I-l.M 

AA(I,J)-A(I,J) 

510  CONTINUE 
500  CONTINUE 
C  WRITE(*,*) 

C  DO  520  I-l.M 

C  WRITE(*,  (1X,<M>F6.2,3X,F6.2) ‘ )  (A(I,J) ,J-1,M) ,B(I) 

C520  CONTINUE 


C 

C--- 

C 

C 

C 

C 

C540 


550 


560 


C 

C 

C 

C570 

C 

c--- 

C 


590 

580 

C 

C 

C 

C600 

530 


DO  530  K— 2 ,M 

IF(IOPT  .NE.  0)  THEN 

--FIND  OPTIMUM  'TOP'  ROW 


PRINT*, 'A: ' 

DO  540  I-l.M 

WRITE (* ,  (10F8 . 3) ' )  (A(I,J) .J-l.M) ,B(I) 
CONTINUE 
Kl-K-1 

AMAX-ABS(A(K-1,K-1>) 

DO  550  KK-K.M 

IF(ABS^A(KK, K-l) )  .GT.  AMAX)  THEN 

AMAX-ABS(A(KK,K-1)) 

ENDIF 

CONTINUE 

IF(Kl.NE.K-l)  THEN 
DO  560  I-K-l.M 
ADUM-A(K-1, I) 

A(K-1 , I)— A(K1 , I) 

A(K1 , I )-ADUM 
CONTINUE 
BDUM-B(K-l) 

B(K- l)-B(Kl) 

B(K1)— BDUM 
ENDIF 
ENDIF 

PRINT*,  'B: 1 
DO  570  I-l.M 

WRITE!*, ' (10F8 . 3) ' )  (A(I,J) ,J-1,M) ,B(I) 
CONTINUE 


DO  580  I-K,M 

B(I^-B(I)^B(K-1)*(A(I,K-1)/A<K-1,K-1.)) 

D°A(I,J)"a(I.J)-A(I,K-1)*(A(K-1,J)/A(K-1,K-1)) 
CONTINUE 
CONTINUE 
PRINT*, ’C: ' 

^°WRITE (*|  * *( IX , <M>F6 . 2 , 3X ,  F6 . 2 )  *  )  (A(I  ,J)  ,J-1,M)  ,B(I) 
CONTINUE 
CONTINUE 


DO  620  J-I+l.M 

B(I)-B(I)-A(I,J)*B(J) 
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non 


ozu 

610 

C630 

C 

650 

C 

640 


CONTINUE 

B(1)-B(I)/A(1.I) 

CONTINUE 

DO  630  I-l.M 

WRITE(* ,*)  (A(I.J) ,J-1,M) ,B(I) 
CONTINUE 


RMAX-0 . 

DO  640  I-l.M 
S— 0 . 


DO  650  J-l.M 

S-S+AA(I,J)*B(J) 

CONTINUE 

R-ABS(S-BB(I)) 

IF(ABS(R)  .GT.  RMAX)  RMAX-R 
WRITE!* , * (1X.G12.5 , 10X.G12 . 5) ' ) 


B(I).  R 


CONTINUE 

IF(RMAX  .GT.  l.E-10) 

+  WRITE!*, ’(A)’)  '***  WARNING,  EQUATION  SOLVER  MAY  BE', 
+  1  INACCURATE  ***' 


RETURN 

END 
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figure  1  —  Experimental  Set-up 
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Figure  2a  —  RTD  Locations  (or  Testing  Frequency  of  1  Hz 
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Rgure  2b  —  RTD  Locations  for  Testing  Frequency  of  <  Hi 
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Figure  B  — Temperature  Profile  1  Hz,  £  =  ±0.52% 
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Figure  6  —  Temperature  Profile  4  Hz,  c  =  ±CM 
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Figure  8  —  Energy  Density  Rate  Histories,  4  Hz,  £  =  ±  0.45% 
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Figure  11  — Total  Energy  Densities  to  Crack  Initiation 


